clc
clear
close all


%% 参数初始化
TC = 1/(480 * 1000 * 4096); % 时延结果采用的时间最小单位
B_subcar = 30e3; % 子载波带宽 Hz
subcar_spacing_comb = 4; % comb间隔，每4个子载波放置1个SRS
fspacing_SRS = subcar_spacing_comb * B_subcar; % SRS信号的频率间隔
num_subcar_SRS = 816; % SRS的有效子载波数


%% 计算与 abs(ifft(Yf./pilot)) 对应的时刻
% SRS信号的时间间隔
tspacing_SRS = 1/fspacing_SRS; 
% SRS的每个子载波对应的时刻
tt = (0:num_subcar_SRS-1)*(tspacing_SRS/num_subcar_SRS);
% SRS的每个子载波对应的归一化时刻
t_TC = tt/TC;

%% 发射的SRS信号频域复数数据 1×816 的向量
load('..\pilot.mat');
% load("..\example_64Tc.mat");
% Yf = example_64Tc;

%% 获取频域数据的文件夹
folder = '../ant_data';  % 文件夹路径
matFiles = dir(fullfile(folder, '*.mat'));  % 获取所有符合要求的文件信息
time_delay = zeros(800,1);
%% 前400份为单天线的频域数据
for ii = 1:400
    YY = get_data(matFiles,folder,ii);
    TC_final = estimate_T(pilot,YY,t_TC);
%     fprintf("TC final = %.2f\n", TC_final);
    time_delay(ii) = TC_final;
end

%% 后400份为多天线的频域数据
for ii = 401:800
%     ii=randi([401 800],1,1);
    YY4 = get_data(matFiles,folder,ii);
    time_delay_ant = zeros(4,1);
    for antj = 1:4
        YY = YY4(antj,:);
        time_delay_ant(antj) = estimate_T(pilot,YY,t_TC);
    end
%     hhhhhhhhhhhhhhh
    time_delay(ii) = time_delay_merge(time_delay_ant);
end

%% 把 time_delay 写入 answer.txt
fileID = fopen('../data/answer.txt', 'w');
for value = time_delay
    fprintf(fileID, "%.2f,\r\n", value);
end
fclose(fileID);




function YY = get_data(matFiles,folder,ii)
filename = fullfile(matFiles(ii).folder, matFiles(ii).name); % 获取文件名及路径
data = load(filename); % 加载MAT文件中的数据
variable_names = who('-file', filename); % 获取MAT文件中的变量名
variable_name = variable_names{1}; % MAT文件中只有1个变量
YY = data.(variable_name); % 获取MAT文件中的变量值 
end


function TC_final = time_delay_merge(time_delay_ant)
TC_final = mean(time_delay_ant);
end



function TC_final = estimate_T(pilot,YY,t_TC)
HH = YY./pilot;
hh_abs = abs(ifft(HH));
% num_subcar_SRS = length(XX); % 信号长度

% hh_abs(locs)=peaks
[peaks, locs] = findpeaks(hh_abs);
% TC_final = t_TC(locs);
% fprintf("TC final = %.2f\n", TC_final);
% 最大峰值
max_peak = max(peaks);
% 筛选峰值：幅度 >= 最大峰值 / K & 位置不在开头
locs_sele = locs(peaks >= max_peak / 3 & locs ~= 1);
locs_sele_first = locs_sele(1);

% 获取峰值及其相邻位置的数据共3个
peak_data = hh_abs(locs_sele_first-1:locs_sele_first+1);
time_data = t_TC(locs_sele_first-1:locs_sele_first+1);
% 使用二次多项式对这3个数据进行拟合
p = polyfit(time_data, peak_data, 2);
% 计算二次多项式的顶点（即真实峰值对应的时刻）
TC_final = - p(2) / (2 * p(1));


% 绘制时域
% figure;
% plot(t_TC, hh_abs);
% grid on;
% xlim([0 256])
% xline(TC_final)
% xlabel('TC')
% ylabel('abs(ifft(Yf./pilot))')
% title(['TC final = ',num2str(TC_final),'TC'])
% hhhhhhhhhhhh
end




